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Abstract 


Me consider the shape of a gas bubble which rises through a quiescent incompressible, Newtonian fluid at 
!nter"-e.^iate Reynolds numbers. Exact numerical solutions for the velocity and pressure fields, as well as the 
bubble shape, are obtained using finite difference techniques and a numerically generated transformation to an 
orthogonal, boundary-fitted coordinate system. No restriction is placed on the allowable magnitude of defor- 
mation. 


Introduction 


In spite of intensive investigation for more than 70 years, the theoretical problem of bubble and drop 
motion in an unbounded viscous Newtonian liqu'd, which may either be quiescent or undergoing some prescribed 
motion far from the bubble or drop, has remaiiied essentially unsolved. The main difficulty, in additio.i to 
the usual nonlinearity of the equations of motion at finite Ryenolds number, is that the bubble or drop shape 
is unknownand required as part of the solution of the problem. As a consequence, the boundary conditions at 
tne bubble or drop surface are nonlinear, and, in addition, the solution depends on the prior history of the 
bubble or drop motion and interface shapes, even in the creeping motion limit. 

This problem of bubble or drop motion i" a viscous liquid, with an unknown boundary shape, is an example 
of the general class of so-called "free b . y" problems of fluid mechanics. Although a number of methods 
exist for this type of problem, which can . ed to analyze the motion of bubbles or drops, all but purely 
numerical methods inevitably s. ffer from sone restriction. Included is the asymptotic technique of "domain 
perturbations" which has beer, applied, fo.‘ example, for buoyancy-driven motions of a drop at finite Reynolds,' 
and for a drop in a general linear "shear" flow,^ but is restricted to small dei rmations from a known (or 
guessed) boundary shape. Boundary integral techniques art .lot restricted in the degree of deformation (and 
thus provide a powerful tool to study bubble and drop deformation'"^), but are only applicable in the limits 
of either creeping flow, or potential flow, where the governing differential equations are linear. For more 
general conditions, this leaves us with numerical methods which are not limited, in principle, either by the 
allowable degree of deformation or by lir arity of the governing equations. Such methods have not been 
applied directly to the problem of calculating bubble or drop shapes in viscous flow so far as wr are aware. 
However, in most other applications to free boundary probletit in fluid mechanics, the numerical r.chods of 
choice have been ased u..on finite element formulations. At least in part, this has been a consequence of the 

loss of accuracy which occurs when finite difference techniques are applied in domains with boundaries that 

are not coincident with coordinate lines or surfaces. Thus, if one considers only the classical orthogonal 
coordinates, such as cylindrical, spherical, etc., the use of finite-difference methods is generally unaccept- 
able for free boundary problems. The present paper explores the alternative possibility of finite-difference 
solutions oased on a numerical method of constructing a system of orthogonal, boundary-fitted coordinates, for 
the problem of streaming flow past a bubble. We do not claim or intend to imply "superiority" in any sense 
over other possible numerical approaches to the same problem. Indeed, the methods described here are not at a 
sufficiently advanced stage of development for such comparisons to be meaningful, even if one were philosophi- 
cally inclined to make them! 

A detailed description of the methods of orthogonal mapping will soon appear in the Journal oj' Computa~ 
tional Physios,^ and a more detailed description of methods and results for the application to the motion of a 

bubble or drop in a viscous fluid is presently in preparation. Here, we simply outline the method of solution 

and present two cxamnles of the solution for streaming flow past a bubble at finite Reynolds number as an 
illustrition of i is application. 


Orthogonal Happing 

The idea which we pursue is thus to obtain orthogonal boundary-f 1 t'.ed coordinates for the domain exterior 
to a bubble whose sh< e is unknown, though smooth and generally nonspheri cal . In the present development, 
the c. . .nown shape is generated via an iterative procedure starting from some initial guess. At each step, 
with the boundary shape specified, mappit^g functions for the boundary-fitted coordinates are generated numer- 
ically, the equations of motion are then solved in the transform (....main and the normal stress balance at the 
bubble surface is used to generate an improved shape. We restrict ourselves to steady, axisyimetric 
configurations and discuss the mapping problem in 20, with the axisymmetric boundary shape generated by rota- 
tion about the axis of synmetry which is thus required to be a coordinate line for the transform coordinates. 

In the remainder of this section, we outline a method of obtaining the desired coordinate transformation. 
From a mathematical point of view, we require a pair of functions x(r‘,n) and y(C,n) wtiich map poi ts of the 
physical domain onto a unit square, 0 S C.t! S 1, in the *ransform dom.ip, with lines of constant C and q 
being orthogonal. For conven i eiice , we designate the bubble surface as ^ » 1 , and the upstream and downstream 
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axes of synmetry as t “ 1 and n “ 0, respectively, with ^ 0 corresponding to infinity. There has, of course, 

been a great deal of recent research aimed at the problem of obtaining numerically generated coordinate map- 
oings. Included in this work are methods based on. the solution of a pair of elliptic equations for the mapping 
functions,’ conformal mapping,* direct integration of "Cauchy-Rlemann"-type equations as an Initial value 
problem starting from a boundary*’** and other methods of orthogonal mapping which are equivalent to conformal 
mapping with less restrictive constraints on the ratio of the diagonal components of the metric tensor (the 
latter are, in fact, most closely related to the present approach) Limitations of space prevent a de- 

tailed review of this prior work. However, in general, the resulting coordinate systems are either nonorth- 
ogonal,’ ill-conditioned in the sense of extreme sensitivity to boundary shape and/or (the possibility of) 
highly nonuniform spacing of coordinate lines (conformal mapping,* some types of "orthogonal mapping'"*!*^) or 
only suitable i some local subdomain (integrations of Cauchy-Riemann equations* • *° ) . 

The present objective is a numerically generated mapping which is applicable in the whole domain, automat- 
ically orthogonal and free of the usual sensitivity problems of conformal mapping. Our basis is conventional 
tensor analysis, yielding equations for x(^,n) and y(C>n) which are coordinate invariant. These equations 
follow almost trivally from the observation that a Cartesian coordinate x (or y) is a linear scalar function 
of position, so that grad x (or grad y) is constant and 

di V grad x ■ 0 (l) 


The latter is nothing more than the covariant Laplace equation for x. When expressed in terms of the desired 
(but as yet unspecified) coordinates, it becomes 


g'-'x. 


•J 


(2) 


in which g'-* is the ij component of the metric tensor and denotes regular covariant differentiation. Al- 
though the solutions of Eq. (1) (and the similar equation for y) will not generally yield orthogonal coordi- 
nates, we further specify that 


q. . » C 

Vir i^J 

(3) 
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with the being "scale" factors for the system. In this case, the equations governing the transfor- 

mation (mapping) functions become 
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and the solution of these equations, subject to appropriate boundary conditions (which we shall discuss in the 
next sectio’’', wi i I yield orthogcnal coordinates for any f, which can thus be chosen freely. It may be noted 
that confer. mapping corresponds to the restrictive choice f(ii,r|) » 1, for all 

In general, the "most appropriate" choice of f depends on the type of mapping required The problem of 
direct mapping with fixed boundary shape and a specified distribution of coordinate nodes along the boundary 
is discussed elsewhere.* Here, we consider only the mapping problem in which the boundary shape is unknown 
and required as part of the solution of the overall problem. In this case, f(5,ri) can be specified directly 
as a function of ^,q, with the form for f chosen so as to yield desired properties of the transform coordi- 
nates (e.g. nonuniform spacing of coordinate lines in some region of the domain). 

The fluid dynamics problem — basic formulation 

Let us now return to the problem of uniform streaming flow past a bubble. In this case, we adopt the very 
simple form for f, f(5,n) ^ In addition, we introduce a relatively simple modification of the mapping 

procedure outlined above to take care of th. lact that infinite values of the mapping functions x and y, 
corresponding directly to an infinite domain, cannot be generated numerically. To avoid this difficulty, we 
simply calculate the mapping fr^.m the unit square in the Ei,ri plane to an auxiliary finite domain, which is 
then transformed to the physical domain by a conformal inversion. 

Now, one great advantage of orthogonal coordinates, in addition to avoiding inaccuracy of numerical 
approximation in ncnorthogonal coordinates. Is that physical components of vectors and tensors can be used 
instead of covariant or contravariant ones. The governing Navier-Stokes equations, plus boundary conditions, 
can thus be expressed in a straightforward manner in terms of the resulting boundary-fitted coordinates 
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obtained by rotation of the two-dimensional c^^.^dinates given by x(^,n) and o(C>n) (where x is parallel 
to the axis of symmetry and a is the distance to this axis along a normal through the point of interest). If 
we introduce the streamfuncrion i|j, and use standard expressions for the invariant differential operators in 
general orthogonal curvilinear coordinates, the Navier-Stokes equations are 


I 3iji 3 (C\ 

h^h^ 3n h^h^ 3C 3n 'O' 

C - 0 

where C is the vorticity, Re ” -r — 


, d is the equivalent diameter of the bubble and 



The streamf unction at infinity, for a . liform streaming flow, takes the form 


Thus, to avoid dealing with large (or infinite) numbers, we actually solve for 
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where 4' is the potential flow solution for flow past a spherical bubble with the given form at infinity, 
i .e. 


^>3 - - 53 ) 


(12) 


Now, Eqs. (4) and (5), rewritten in terras of I' , are to be solved for 4' , Z and the bubble shape subject 
to the boundary conditions 

is bounded, 4 =• 0; at infinity (i.e. 5*0) 

ij « 0, 4 = 0; at n = 0, 0 = 1 (symmetry axis) 

and, at the bubble surface, 
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(14) 

(zero normal velocity) 1 

(15) 

(zero tang, stress) \ at 5 “ 1 

(16) 

(normal stress balance) / 

(17) 


The first term in (16) is the hydrostatic pressure; is the drag coefficient; p^ is the dynamic pressure 
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dPv: 


u, is the surface velocity; i is the normal component of viscous stress at the surface. We » 


; Y is 
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the surface tension, and '^’(„)^ Rod "^(n) *'"* curvatures in two perpendicular directions. 


Numerical scheme 


In order to solve Eqs. (7) and (8) of the preceding section, together with Eqs. (5a) and (5b) for the 
mapping function , x(5,h) and o(5,n), we used a uniform 41x41 grid in the domain, 0 < 5,n < 1. The computa- 
tions were carried out using single precision arithmetic on a VAX-II computing system, which has a round-off 

error of 0(10*8). Tnus, with an 0{h^) finite-difference scheme, this mesh size represents the practical 
limits of resolution in order ti" r truncation error be comparable to this rounu off error divided by h^ (when 
computing second derivatives). 

The numerical scheme itself must be fast, highly stable and applicable to elliptic equations of quite 
general form. In the work reported here, we adopt the ADI scheme of Peaceman and Rachfcrd and treat all 

equations of the problem (i.e. the equations of motion for 4 end <i', and the two mapping equations for x and 

o) as "quasi -time-dependent", by writing them in the standard form 
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with 3 / 3 t representing a "fictitious" (or artificial) time derivative as required by ADI. 
the iteration parameter (i.e. time step) was determined^ '* to be 0(h). 


An optimal value of 


Boundary conditions for Eqs. (3), (5) and (6) are straightforward [see Eqs. (13)-06) plus Ref. 6], with 
the exception fo conditions at the bubble surface. Here, the necessary boundary values of vorticity are cal- 
culated indirectly from the boundary condition (16) on tangential stress usinq a natural extension of the 
method for a solid boundary suggested by Dorodnitsyn and Heller*^ and Israeli^® and utilized previously for a 
spherical drop.'^ At each new iteration, say n, the new value of the boundary vorticity r." is determined from 
its previous value and the previous value of the tanqential stress, as 
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where the optimal b was found (by trial and error) to be approximately 0.2. When the solution has converged, 
of course, the tangential stress will be zero. Boundary conditions for x(C,r;) and o(f.,ii) at 5 • I must also 
be discussed briefly. Both x and a cannot be specified directly at ■ 1 if the condition “ 0 is satis- 
fied (i.e. the coordinates are to be orthogonal) as the problems for x and u are then overdetermined. We 
would, on the other hand, like to approach the final solution for bubble shape iteratively starting from some 
initial guess. This involves incrementing the bubble boundary to create a new shape at each Iteration, based 
upon the normal stress imbalance at the interface at the preceding iteration. However, in view of the re- 
striction on simultaneous specification of x and the necessary small displacement of the bubble boundary 
must be carried out indirectly rathe'" than specifying increments in x(1,n) and ^'(l,n) directly. This is ac- 
complished by changing the mapping Itself (rather than the position of the bubble surface) via incremental 
changes in the scale factor h_. of the mapping, i.e. 
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The incremented h; 


'.«1 


is then used to generate "equivalent" boundary conditions for 
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The normal stress difference, 3, has to be normalized before it is used in (21) for changing the bubble shape 
because of the indeterminacy due to incompressibility (.3 contains an integration constant): this indeterminacy 
is removed by requiring that the volume of the bubble remain con. .ant. 


’’he overall solution algorithm may thus be schematically represented as follows; 

(1) Start with an Initial guess of the shape. Here we choose a spherical shape, i.e. a circle in a plane 
through the axis of symmetry. Hence, with f(3,n) = as indicated earlier, the mapping is initially x - 
3cos”n and y - Tsincii, corresponding to polar coordinates in the plane through the axis of symmetry. 

(2) For the given bubble shape and coordinate mapping, compute a new approximation for the dynamic fields 
(3 and ',) by advancing the solution of the Navier-Stokes equations one iteration (I.e. one ADI time step). 

(3) Calculate the normal stress terms at the bubble surface, and if the condition (17) is not satisfied, 

increment the bubble shape by a small amount by incrementing ir(l,h) using Eq. (21), and obtaining correspon- 
ding boundary conditions for • 

(4) Calculate a new orthogonal mapping fitting the new bubble shape by solving Eqs. (5a, b) with appropri- 
ate boundary conditions (in practice we do only one ADI iteration on the mapping equations). 

(5) Repeat this process starting with step (2) until convergence is achieved. 


Results 


We consider two cases here of streaming flow past a bubble. 

Case A; Rc - 2.47, W- - 4.00 
'•ase B: Re - 100. , We - 2.00 
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The final bubble shape Is depicted for Case A In Fig. t, where we show a portion of the final coordinate mesh 
and the upper half of the bubble boundary In the plane through the axis of symmetry. The flow Is from left to 
right. The corresoonding streamlines and lines of constant vorticity are shown in Figs. 2 and 3. It may be 
noted that the bubble shape Is in qualitative agreement with available experimental results. “ Indeed, the 
drag coefficient calculated here is 9.17i whereas the measured value at the same Reynolds number but somewhat 
larger Ue was 9.37. It may be noted that the drag coefficient was found experimental ly‘* to be Insensitive to 
We for large We. The streamlines and lines of constant vorticity for Case B are shown In Figs. It and 5, from 
which the bubble shape can also be discerned. Again, the flow Is from left to right. It Is thus evident that 
the bubble Is actually flattened to a greater degree in the front and Is more rounded at the rear. Shapes of 
this general type have been previously observed experimentally for similar values of Re and Ue,‘^ although a 
shape which is rounded in the front and flattened at the rear, which will occur for larger Reynolds number or 
larger Weber number, I .e. smaller surface tension. Is much more common. Each example required about an hour 
of CPU time on a VAX-ll computer, starting from the irrotational flow pas a sphere as an initial guess in 
both cases. The cost is thus on the order of $10. 
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Figure 1. Coordinate mesh for Case A. 
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